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Abstract 

With Monte Carlo methods we investigate the dynamic relaxation of the fully 
frustrated XY model in two dimensions below or at the Kosterlitz-Thouless 
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and universality is confirmed. The critical exponent is measured for different 
temperature and with different algorithms. 
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I. INTRODUCTION 



There has been a long history for the study of universal scaling behaviour for critical 
dynamics. It is well known that there exists a universal scaling form when dynamic systems 
almost reach equilibrium or in the long-time regime of the dynamic evolution. When a 
dynamic magnetic system has evolved for a sufficiently long time, the magnetization decays 
exponentially. The characteristic time scale for this regime is t T ~ t~ uz or ti ~ L z with r 
being the reduced temperature and L being the lattice size. Before this exponential decay 
the time evolution of the magnetization obeys a power law t~^l vz . All this universal scaling 
behaviour can be characterized by a set of three critical exponents, two static exponents /3, 
v and one dynamic exponent z. 

Is there universal behaviour in the short-time regime of the dynamic evolution? For long, 
the answer had been no. It was believed that the behaviour of the dynamic systems in the 
short-time regime depends essentially on the microscopic details. However, it has recently 
been discovered that universal scaling behaviour emerges already in the macroscopic short- 
time regime, starting at a time just greater than a microscopic time scale t mic [p]]. Important 
is that extra critical exponents should be introduced to describe the dependence of the scaling 
behaviour on the initial conditions or to characterize the scaling behaviour of special 
dynamic observables A typical dynamic process is that a magnetic system initially 

at high temperature with a small initial magnetization, is suddenly quenched to the critical 
temperature (without any external magnetic field) and then released to a time evolution 
with a dynamics of model A. A new dynamic exponent xq has been introduced to describe 
the scaling dimension of the initial magnetization. More surprisingly, at the beginning of 
the time evolution the magnetization undergoes a critical initial increase |l],|6],[7| 

M(t)~m t e (1) 

where 9 is related to the scaling dimension x of the initial magnetization m by 9 = 
(x - P/u)/z. 

The physical background for the initial increase of the magnetization has not been com- 
pletely clear. The critical initial increase has little to do with the symmetry breaking below 
the critical temperature. Actually similar phenomena can be observed even in case that the 
magnetization is not an order parameter. For the 6-state clock model and the XY model, 
below the Kosterlitz Thouless phase transition temperature no real long range order ap- 
pears. The normal magnetization is not an order parameter. However, the power law initial 
increase is still observed in numerical simulations |||||. 

On the other hand, in recent years much attention has been drawn to statistical systems 
with frustration or quenched randomness. Critical behaviour of these systems is often quite 
different from that of regular systems. Due to the severe critical slowing down, numeri- 
cal simulations of these systems are extremely difficult. Especially, our knowledge on the 
dynamic properties of these systems is poor. 

In this paper, as a first approach to the short-time dynamics of statistical systems with 
frustrations, we investigate numerically the dynamic fully frustrated XY model below and 
at the Kosterlitz-Thouless phase transition temperature. We concentrate our attention on 
the scaling behaviour of the magnetization and its dependence on the initial value. Much 
effort is devoted to the understanding of the special properties induced by the sublattice 



2 



structure of the ground states. An investigation of the universal short-time behaviour of 
the dynamic systems is not only conceptually important, but also in the practical sense. 
Numerical simulations for the Ising and Potts model show that we may obtain already from 
the short-time dynamics the static critical exponents as well as the dynamic exponent z, 
which are normally defined and measured in equilibrium or in the long-time regime of the 
dynamic evolution |7|, 10|, 11], g|,I^,fL3l. Since our measurements in the short-time dynamics 



are always carried out at the beginning of the time evolution and the average is really a 
sample average rather than a time average based on the ergodicity assumption as in the 
measurements in equilibrium, the dynamic approach might be free of critical slowing down. 

In the next section, the fully frustrated XY model is briefly introduced. In Sec. 3 and 
Sec. 4, numerical results are reported for the dynamic fully frustrated XY model. Finally 
we give some conclusions. 

II. THE FULLY FRUSTRATED XY MODEL 

The fully frustrated XY model in two dimensions (FFXY) can be defined by the Hamil- 
tonian 

H = K £ S> ■ S, , (2) 

<ij> 

where Si = (S^x, S^) is a planar unit vector at site % and the sum is over the nearest 
neighbours. In our notation the inverse temperature has been absorbed in the coupling K. 
Here fij takes the values +1 or —1, depending on the links. A simple realization of the 
FFXY model is defined by taking = — 1 on half of the vertical links (negative links) while 
+1 on the others (positive links). This is shown in Fig. [I]. The links marked by dotted lines 
represent the negative links. 

In the FFXY model two phase transitions exist, the Kosterlitz-Thouless phase transition 
(XY-like) and the second order phase transition (Ising-like). This is very different from 
the regular XY model. Much effort has been made to locate the critical points for both 
transitions and measure the corresponding critical exponents. In a recent paper [H]], it is 



reported that the XY-like phase transition temperature is Tkt = I/A'xt = 0.446 while the 
Ising-like phase transition temperature T c = 1/ K c = 0.452. These two measured values for 
the phase transition temperature are slightly different. Earlier results for the FFXY model 
and its ordering dynamics can be found in references []I5-2"0]. 



One of the most important properties of the fully frustrated XY model is the sublattice 
structure of its ground states. In the regular XY model, in the ground state all spins orient 
in the same direction. However, for the FFXY model the lattice should be divided into four 
sublattices. All spins in each sublattice orient in the same direction, but the directions for 
the four sublattices are different and connected in a certain way, as it is shown in Fig. [I]. For 
the whole lattice, the ground state preserves the global 0(2) symmetry at any temperature, 
i.e. spins can rotate globally. No real long range XY-like order emerges in the FFXY model. 
This situation is similar as for the regular XY model. Besides the ground state shown in 
Fig. [TJ, there is another ground state which is just obtained by translating the configuration 
in Fig. |I] by one lattice spacing in y direction. Below the Ising-like critical temperature T c , 
the Z 2 symmetry is broken. A second order phase transition occurs. 
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In this paper we study the short-time dynamic properties of the FFXY model below 
or at the XY-like phase transition temperature Tkt- Below the critical temperature Tkt-, 
the FFXY model remains critical in the sense that the correlation length keeps divergent. 
Therefore, a scaling form is expected even below the critical temperature Tkt- However, 
the critical exponents may vary with respect to the temperature. Such a phenomenon has 
been observed in the 6-state clock model and the regular XY model |||§. The dynamic 
properties related to the Ising-like phase transition will not be discussed in this paper. 

III. SUBLATTICE STRUCTURE OF DYNAMIC EVOLUTION 

Let us consider a dynamic relaxation process starting from an initial state with a very 
high temperature and small magnetization As a direct generalization of the XY model, 
the magnetization for the FFXY model may also be defined as 

^(*) = ^E3 ( 3 ) 

with L being the lattice size. As in the XY model, the magnetization here is also not an order 
parameter. To achieve an initial magnetization, we introduce an initial external magnetic 
field, e.g. in the x direction, as in the numerical simulation of the XY model ||. Taking 
into account that the initial state is at a very high temperature, the initial Hamiltonian can 
be written as 

H = 2h^S iiX . (4) 

i 

To prepare the initial state, we update the system described by the initial Hamiltonian H 
until it reaches equilibrium. Then the generated configurations of this initial system are 
used as initial configurations of the dynamic system. The initial magnetization generated is 

M(0) = (mo, 0) ss (h, 0), h -> 0. (5) 

There are, of course, many other methods to generate an initial magnetization. However, 
it has been demonstrated that the universal behaviour does not depend on these microscopic 
details, i.e. how a magnetization is constructed ||. Effects of the microscopic details of the 
initial configurations are swept away in almost one Monte Carlo time step. 

After preparation of an initial configuration, the system is released to the dynamic evolu- 
tion of model A below or at the XY-like transition temperature. In this paper, the Metropolis 
algorithm is mainly used. To confirm universality, however, some simulations are repeated 
with the heat-bath algorithm. We stop updating the dynamic system at 150 Monte Carlo 
time steps, and repeat the procedure with another initial configuration. The average is over 
the independent initial configurations as well as the random numbers. We have performed 
simulations with lattice sizes L = 8, 16, 32, 64 and 128. The total number of samples 
for the average is between 30 000 and 60 000, depending on lattice sizes, initial states and 
algorithms. Errors are estimated by dividing the samples into three or four groups. 

In Fig. 0, the time evolution of the magnetization is displayed in double-log scale with a 
solid line for the lattice size L = 64 and the initial magnetization m = 0.02. The temper- 
ature is taken to be T = 0.400, which is slightly below the XY-like transition temperature 
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T KT = 0.446 given in reference |L4f (or T KT = 0.440(2) in p|). In the figure M(t) is 



the x-component of the magnetization M(t). The ^-component of the magnetization M(t) 
remains zero since the initial value is zero. From the figure we see that the magnetization in- 
deed increases after some time steps. The universal power law behaviour becomes apparent 
after about 20 — 30 Monte Carlo time steps. From the slope of the curve, one measures the 
exponent 6 = 0.184(6). In this simulation a global uniform initial external magnetic field h 
has been applied to the whole lattice, or in other words, the initial magnetization density 
distribution is uniformly generated. We call this initial state the global start. 

Does the sublattice structure of the ground state play some role in the dynamic evolution? 
In equilibrium, it is known that the XY susceptibility should be calculated separately for 
each sublattice since the orientations of the spins in different sublattices differ from each 
other. But further understanding of the sublattice structure in numerical simulations can 
not so easily be achieved due to the 0(2) symmetry and large fluctuations. In the short- 
time dynamics, the situation is somewhat different. The 0(2) symmetry is violated by 
the initial magnetization. One can really measure the time evolution of the magnetization 
for each sublattice shown in Fig. |I] separately. The results are also included in Fig. 
The two upper dotted lines represent the time evolution of the magnetization for the two 
sublattices connected to the positive links while the lower two dotted lines are those for the 
two sublattices connected to the negative links. In is very interesting that for sublattices 
connected to the positive links the magnetization increases from the beginning while for 
sublattices connected to the negative links they drop significantly in the first time steps. 
However, after a certain number of time steps the magnetization tends in all cases to the 
same universal power law behaviour even though the magnitudes remain different for the 
sublattices on the positive links and the negative links. 

What is the relation between the magnetizations on the positive links and the negative 
links? In Fig. |3], the ratio r(t) of the magnetization averaged over two sublattices on the 
positive links and that on the negative links is plotted with a solid line. This ratio stabilizes 
to a constant very quickly within 10 time steps. Averaging this ratio in a time interval 
[10, 150], we get r = 2.4147(8). What is this ratio? In the ground state drawn in Fig. [I], the 
magnitude of the angles between spins on the positive links and the x axis is tt/8, while that 
between spins on the negative links and the x axis is 37r/8. The ratio r should be nothing but 
the ratio of the x-components of the spins on both the positive links and the negative links. 
We see this by calculating r th = cos(-7r/8)/ cos(37r/8) ~ 2.4142. The consistence between r th 
and the measured one r = 2.4147(8) is remarkable. 

For a better understanding of this point we now divide the lattice into two sublattices, 
the sublattice on the positive links and the one on the negative links. We have performed 
a simulation with different initial magnetizations for the two sublattices, m 0p and m 0n re- 
spectively, taking, for example, m 0p /m 0n = r th and keeping the global initial magnetization 
m = 0.02. We call such an initial state two sublattice start. The time-dependent mag- 
netizations with this initial condition are plotted in Fig. |j. The solid line is the global 
magnetization while the upper and lower dashed line represent those for sublattices on the 
positive links and negative links respectively. Now, the decline for the magnetization on the 
negative links has disappeared. The corresponding ratio r(t) is plotted in Fig. |3| with the 
dotted line. It is a constant from almost the very beginning of the time evolution. The 
averaged value is r = 2.4146(6). However, the time period for the magnetization to enter 
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the universal power law behaviour is again 20 — 30, almost the same as that in Fig. |2|. 

Before going further we would like to mention the problem of the sharp preparation of 
the initial magnetization. If the lattice size is infinite, in each initial configuration generated 
by the initial Hamiltonian H an exact value (m ,0) of the initial magnetization M(0) is 
automatically achieved. However, in practice the lattice size is finite and the initial magne- 
tization M(0) fluctuates around (mo, 0). This is a kind of extra finite size effect. It causes a 
problem in high precision measurements. In order to reduce this effect, a sharp preparation 
technique has been introduced |]6|,[7|j21[1 . How important the sharp preparation technique is 
depends on the initial magnetization m and what kind of observables one measures. The 
smaller the initial magnetization tuq is, the more important the sharp preparation technique 
becomes. In Fig. f|, the result for the time evolution of the magnetization with the sharp 
preparation technique has been plotted with the dotted line. The curve almost overlaps with 
that without the sharp preparation technique. This shows that the extra finite size effect 
here is already quite small for the lattice size L = 64. For simplicity, our simulations are 
always carried out with no sharp preparation of the initial magnetization. 

Now we go a step further. From the numerical simulations shown in Fig. |2|, [3] and |], and 
the discussions above, we understood that the universal behaviour of the dynamic system 
has close relation with the structure of the ground states. If the initial state is a completely 
random state with zero initial magnetization, the probabilities for the magnetization to 
evolve to different directions are the same, the averaged magnetization remains zero. How- 
ever, if a non-zero initial magnetization is given to a certain direction, the time dependent 
magnetization grows in this direction since the energy is in favour of it. To clarify this point, 
we start with an initial state which is even closer to the ground state. We give different 
orientations to the initial magnetizations for the four sublattices as shown in Fig. [I]. All 
the magnitudes of the initial magnetizations are m = 0.02. This initial state we call four 
sublattice start. The time evolution of the magnetizations for different sublattices are shown 
in Fig. |5|. Here M(t) is just the projection of the magnetization on the direction of the initial 
magnetization. The solid line and dotted line are the magnetizations for the sublattices on 
the positive links while the circles and crosses are those on the negative links. They collapse 
on a line and show completely the same universal behaviour. 

In Table |, values for the critical exponent 9 measured in a time interval [40, 150] for the 
three different initial states are collected. Within the statistical errors they are consistent. 
In the next section we will come back to this point. 

IV. UNIVERSALITY AND SCALING 

In the numerical measurements of the critical exponent 9 we should pay attention to 
two possible effects, the finite size effect and the finite mo effect. For the dynamic process 
discussed in the last section, there are two kinds of finite size effects. One is the extra finite 
size effect from the initial configurations. This extra finite size effect is closely related to the 
problem of the sharp preparation of the initial configurations and has been discussed in the 
last section. 

Another kind of finite size effect is the normal finite size effect which takes place in a 
time scale t^ ~ L z . Whenever the system evolves into this time regime, the magnetization 
will decay by an exponential law exp(— t/ti). In order to see the normal finite size effect, 
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we have plotted in double-log scale the time evolution of the magnetization for the lattice 
sizes L = 8, 16, 32, 64 and 128 with initial magnetization m = 0.02 in Fig. ||]. The upper 
solid line is the time-dependent magnetization for L = 64, while the dotted lines are those 
for L = 8, 16, 32 and 128 respectively. The curves for L = 64 and L = 128 more or less 
overlap. The values for the critical exponent 9 are 9 = 0.181(5) and 0.182(3) for L = 64 and 
L = 128 respectively. Therefore we conclude the finite size effect here is already negligibly 
small for the lattice size L = 64. 

Rigorously speaking, the critical exponent is defined in the limit m = 0. However, it 
is practically only possible to perform the measurement with finite m . The exponent 9 
measured may show some dependence on the initial magnetization rriQ. In general, a linear 
extrapolation to the limit of m = should be carried out 0)0- For this purpose, we 



have performed another simulation for the lattice size L = 64 with an initial magnetization 
m = 0.01. The time- dependent magnetization is also displayed with the lower solid line in 
Fig. The exponent 9 obtained for m = 0.02 and m = 0.01 are 9 = 0.181(5) and 0.179(7) 
respectively. Within the statistical errors they cover each other. Since the dependence of 
the exponent 9 on tuq is already rather weak here, a linear extrapolation is not necessary. 
This is also the reason why the results with and without the sharp preparation of the initial 
magnetization are not so different. In the further discussions in this paper we will work 
mostly with L = 64 and m = 0.02. 

Is the power law scaling behaviour in Eq. (|l]) for the magnetization really universal? 
For example, can it depend on the microscopic details of the initial state, the algorithms 
and the lattice types and even the additional non-nearest interactions and so on? Many 



discussions of this kind have recently been made pl|j22]j9|I41,p3| . Here we have also repeated 



some calculations with the heat-bath algorithm to confirm universality. The result for L = 64 
and m = 0.01 with the two sublattice start is plotted in Fig. [7| with the solid line. The 
initial magnetization mo p and m§ n are given to the sublattices on the positive links and the 
negative links respectively, with m p/wo n = r t h = 2.4142. The long dashed, dashed and 
dotted line are the results with the Metropolis algorithm with itlq = 0.02 and the global 
start, two sublattice start and the four sublattice start respectively. All the measured values 
for the corresponding exponent 9 are listed in Table |. The fact that all the four curves give 
consistent results for the exponent 9 within the statistical errors provides strong support for 
universality. Especially the independence of the exponent 9 on the initial states indicates 
that some physical mechanism closely related to the ground states essentially governs the 
time evolution of the dynamic systems. 

To understand universality in the short-time dynamics, we keep in mind that there are 
two very different time scales in the dynamic systems, the microscopic time scale and the 
macroscopic time scale. A typical macroscopic time scale is t T or t^. Universal behaviour 
emerges only after a sufficiently long time period in the microscopic sense. The time period 
which the dynamic system needs to sweep away the microscopic behaviour is called t mic . One 
expects that the time scale t m i C is still very short in the macroscopic sense. In numerical 
simulations, a Monte Carlo time step can be considered as a typical microscopic unit. Most 
of the numerical simulations for dynamic systems show that t m i C ~ 10 — 50. In Fig. ^ we 
see that this is also the case for the FFXY model. These results are reasonable. Compared 
with the typical macroscopic time scale t T ot £l, the microscopic time scale t mic observed 
in numerical simulations is indeed very small. In some cases, universal behaviour emerges 
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actually already in one Monte Carlo time step, e.g. in the numerical measurement of the 
critical exponent 9 for the two-dimensional Potts model [0,^3]- Such a clean behaviour in 
the very short-time regime is somehow unexpected. 

It is well known that for the XY-like phase transition no real long order emerges even 
below the transition temperature Tkt- The system remains critical in the sense that the 
correlation length is divergent. A similar scaling behaviour is expected for any temperature 
below Tkt- However, the critical exponents may vary with respect to the temperature. In 
Fig. H the time evolution of the magnetization for different temperatures has been displayed. 
We clearly see that there exists power law behaviour for all the temperatures below Tkt- 
However, the critical exponent 9 varies essentially as can be seen from the results given 
in Table ||. The dependence of the critical exponent 9 on the temperature in the FFXY 
model is qualitatively the same but stronger than that for the regular XY model M. As the 
temperature decreases, the exponent 9 first increases and then slowly decreases. Around 
the XY-like transition temperature Tkt — 0.446 |TJJ] (or 0.440(2) in reference [|HJ), the 
exponent 9 becomes quite small. Therefore, the power law behaviour around the critical 
temperature Tkt is less prominent compared with that at the lower temperature. We have 
also performed some simulations with temperatures above Tkt- Since the correlation length 
is divergent exponentially when the temperature approaches the XY-like transition point 
from above, however, no rigorous information for the XY-like transition temperature Tkt 
could be obtained from our data. Further investigation along this direction is needed. 

V. CONCLUSIONS 

We have numerically simulated the dynamic relaxation process of the fully frustrated 
XY model (FFXY) in two dimensions starting from an initial state with a very high tem- 
perature and small initial magnetization. Special attention has been drawn to the sublattice 
structure of dynamic evolution induced by the ground state. Universal power law behaviour 
is found independently of the sublattice structure of the initial states. The critical exponent 
9 has been measured for different temperatures below Tkt- Universality has been further 
confirmed by carrying out the simulations with both the Metropolis and the heat-bath al- 
gorithm. Many important problems remain open in connection with the present paper, e.g. 
the short-time behaviour of the dynamic FFXY model with respect to the Ising-like phase 
transition or the determination of the critical point as well as other critical exponents from 
the short-time dynamics. 
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TABLES 





Metropolis 


Heatbath 




I 


II 


III 


II 


e 


.184(6) 


.182(5) 


.181(5) 


.186(6) 



TABLE I. The exponent 9 measured for lattice size L = 64 with different types of initial 
configurations and algorithms. I, II and III represent initial states of global start, two sublattice 
start and four sublattice start. 



T 


0.446 


0.440 


0.420 


0.400 


0.350 


0.300 


0.250 


9 


.060(7) . 079(4) .141(5) 


.181(5) 


.245(3) 


.263(2) .260(2) 



TABLE II. The exponent 9 measured for different temperatures with the Metropolis algorithm. 
The lattice size is L = 64. 
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FIGURES 




FIG. 1. One of the ground states for the FFXY model. 
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FIG. 2. The time evolution of the magnetization with the global start for L = 64 and 
rriQ = 0.02 with the Metropolis algorithm is plotted in doubledog scale. M{t) is the x-component 
of the magnetization M(t). The solid line represents the global magnetization while dotted lines 
are those for the four sublattices. 
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FIG. 3. The ratio r(t) for different initial states. 
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FIG. 4. The time evolution of the magnetization with the two sublattice start for L = 64 and 
rriQ = 0.02 with the Metropolis algorithm is plotted in double-log scale. M{t) is the x-component 
of the magnetization M(t). The solid line represents the global magnetization while dashed lines 
are those for the sublattices. For comparison, the dotted line is the global magnetization with the 
sharp preparation of the initial magnetizations. 
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FIG. 5. The time evolution of the magnetization with the four sublattice start for L = 64 and 
mo = 0.02 with the Metropolis algorithm is plotted in double-log scale. M{t) is the projection of 
the magnetization M(t) on the initial direction. The solid and dotted line are the magnetizations 
for the sublattices on the positive links while the circles and crosses represent those on the negative 
links. 
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FIG. 6. The time evolution of the magnetization with the four sublattice start for different 
L and tuq with the Metropolis algorithm is plotted in double-log scale. M{t) is the projection of 
the magnetization M(t) on the initial direction. The solid line represents the magnetization for 
L = 64 while the dotted lines are those for other lattice sizes. 
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FIG. 7. The time evolution of the magnetization for different initial states and algorithms 
is plotted in double-log scale. M{t) is the projection of the magnetization M{t) on the initial 
direction. The solid line represents the mag- netization for L = 64 and ttiq = 0.01 with the 
heat-bath algorithm and the two sublattice start while long dashed, dashed and dotted lines are 
those for L = 64 and mo = 0.02 with the Metropolis algorithm for the global, two sublattice and 
four sublattice start respectively. 
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FIG. 8. The time evolution of the magnetization for L = 64 and tuq = 0.02 with different 
temperatures is plotted in double-log scale. The four sublattice start and the Metropolis algorithm 
are used in the simulations. M(t) is the projection of the magnetization M(t) on the initial 
direction. From above, the temperature is T = 0.250, 0.300, 0.350, 0.400, 0.420, 0.440 and 0.446 
respectively. 
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